## Congressional Bargaining and the Distribution of Grants
## Replication Code for Figures
## Leah Rosenstiel, November 2022

library(tidyverse)
library(EnvStats)
rm(list=ls())

### Read in data ####
floor_amendments <- read_csv("floor_amendments.csv")
floor_votes <- read_csv("floor_votes.csv")
state_pop <- read_csv("state_population.csv")


### Figure 5: Distribution of Winning Coalition Size ####
floor_votes_temp <- floor_votes %>%
  group_by(proposal_id, passed) %>%
  summarise(wc_size = sum(vote,na.rm=T),
            wc_type = "Roll Call Votes")
floor_amendments %>%
  filter(!is.na(share_diff) & !is.na(diff)) %>%
  group_by(proposal_id, passed) %>%
  summarise(wc_size_share = sum(share_diff>=0),
            wc_size_dollar = sum(diff>=0)) %>%
  pivot_longer(cols=starts_with("wc_size"), names_prefix = "wc_size_",
               names_to = "wc_type", values_to = "wc_size") %>%
  filter(wc_type!="dollar") %>%
  mutate(wc_type="Grant Share Increase") %>%
  bind_rows(floor_votes_temp) %>%
  ggplot(aes(x=wc_size, colour=factor(passed), fill=factor(passed))) +
  facet_grid(cols = vars(wc_type),scales = "free_x") +
  geom_density(alpha=0.5) +
  xlab("\nWinning Coalition Size") +
  ylab("Density\n")+
  scale_color_discrete(labels=c("Did Not Passed     ", "Passed"), 
                       name="Amendment Outcome:     ") +
  scale_fill_discrete(labels=c("Did Not Passed     ", "Passed"), 
                      name="Amendment Outcome:     ") +
  theme_minimal() +
  theme(legend.position = "bottom",
        text = element_text(size=10))

### Figure A1: State Population, 2018 ####
state_pop %>%
  filter(year==2018) %>%
  ggplot(aes(x=pop/1000000)) +
  geom_histogram() +
  ylab("Count of States") +
  scale_x_continuous(name="State Population in 2018 (in Millions)") +
  theme_minimal()

### Figure A2: Increases in Population Across States, 1970 to 2018 ####
state_pop_2018 <- filter(state_pop, year==2018) %>%
  select(pop_2018 = pop, state)
state_pop %>%
  left_join(state_pop_2018) %>%
  mutate(pop_increase = (pop_2018 - pop)/pop) %>%
  filter(year==1970)%>%
  ggplot(aes(x=pop_increase)) +
  geom_histogram() +
  scale_y_continuous(name="Count of States",breaks=c(0,5,10)) +
  scale_x_continuous(name="Population Percentage Increase, 1970 to 2018",
                     labels=scales::percent) +
  theme_minimal()


### Figure A3: Amendments to Formula Grant Programs Over Time ####
proposal_level <- floor_amendments %>%
  select(congress, proposal_id, passed) %>%
  unique %>%
  group_by(congress) %>%
  summarize(count_proposals = n(),
         count_passed = sum(passed),
         count_failed = sum(1-passed)) %>%
  ungroup 
ggplot(proposal_level) +
  geom_point(aes(x=congress,y=count_passed,colour="Passed")) +
  geom_line(aes(x=congress,y=count_passed,colour="Passed")) +
  geom_point(aes(x=congress,y=count_failed,colour="Did Not Pass")) +
  geom_line(aes(x=congress,y=count_failed,colour="Did Not Pass")) +
  geom_text(data = data.frame(x=114, y=c(1,2), 
                              text=c("Passed", "Did Not Pass")),
            aes(x=x,y=y,label=text), inherit.aes = F, hjust=-0.1) +
  scale_y_continuous(name="Count of Amendments \n",limits = c(0,20)) +
  scale_x_continuous(name="\n Congress",limits = c(80,118)) +
  scale_color_manual(values=c("lightblue","darkblue"),guide="none") +
  theme_minimal()


### Figure A4: Simulated Reversion to the Mean ####

set.seed(2020)
sim_data <- floor_amendments %>%
  mutate(share_cl = share_cl / 100,
         share_new = share_new / 100,
         share_diff = share_new - share_cl) %>%
  group_by(state_congress) %>%
  mutate(sq_mean = mean(share_cl,na.rm=T)) %>%
  ungroup %>%
  select(share_cl, sq_mean, share_new, state_congress)
sim_data$winner <- ifelse(sim_data$share_new > sim_data$share_cl,1,0)
actual_model <- felm(winner~share_cl|state_congress,data=sim_data)
sim_coefs <- rep(NA,1000)
sim_grants <- rep(NA,1000)
mean_grant_new <- mean(sim_data$share_new, na.rm=T)
var_grant_new <- var(sim_data$share_new, na.rm=T)
beta_param <- EnvStats::ebeta(sim_data$share_new,method="mme")$parameters
for (i in 1:1000) {
  temp_data <- sim_data
  temp_data$share_new_sim <- rbeta(nrow(sim_data), beta_param[1], beta_param[2])
  temp_data$winner_sim <- ifelse(temp_data$share_new_sim > temp_data$share_cl,1,0)
  sim_model <- lfe::felm(winner_sim~share_cl|state_congress,data=temp_data)
  sim_coefs[i] <- sim_model$coefficients[1]
}
as.data.frame(sim_coefs) %>%
  ggplot(aes(x=sim_coefs)) +
  geom_histogram() +
  xlab("Coefficients from Simulation") +
  ylab("Count") +
  geom_vline(xintercept = actual_model$coefficients[1], colour="red") +
  annotate("text", x=actual_model$coefficients[1],y=85, label="Observed = -4.054",
           hjust=-.05,colour="red") +
  theme_minimal()

# How many are more extreme than observed?
table(sim_coefs>=actual_model$coefficients[1])


